Phân tích hồi quy tuyến tính

Premium service Source Web

1 Định nghĩa về hồi quy tuyến tính

Hồi quy (regression) là cách tiếp cận xây dựng mô hình toán học (công thức toán) để biểu diễn mối quan hệ giữa một biến \(y\) với một hay nhiều biến \(x\). Theo nghĩa rộng thì đặc tính của các biến \(x\) và \(y\) có thể là biến liên tục (continuous variable), biến số đếm (discrete quantitative variable) hay biến phân loại (discrete qualitative hay categorical varaible).

Hồi quy tuyến tính đơn biến là dạng căn bản và đơn giản nhất của các mô hình hồi quy.

Hồi quy tuyến tính (simple linear regression model - SLR) theo cách hiểu thông thường \(y = ax + b\) hay \(y = b_0 + b_1x\) là mô hình toán học biểu diễn mối quan hệ giữa hai biến định lượng (continuous variable) là \(x\) và \(y\).

Trong đó biến \(x\) gọi là biến độc lập (independent), biến tiên lượng (predictor), biến giải thích (explanatory) hay là biến hồi quy (regressor), còn biến \(y\) là biến phụ thuộc (dependent), biến dự báo/dự đoán (predicted), biến phản hồi (response) hay là biến hồi quy phụ thuộc (regressand).

Mô hình hồi quy tuyến tính này, có tên gọi cụ thể là hồi quy tuyến tính đơn biến đa thức bậc 1 khi chỉ sử dụng 1 biến \(x\) (bậc 1, lũy thừa 1) để dự đoán cho biến \(y\). Tên gọi mô hình này là Ordinary least squares (OLS) hay là hồi quy tuyến tính theo phương pháp bình phương tối thiểu theo kiểu bình thường (ordinary) là một loại trong các mô hình hồi quy tuyến tính linear regression¹.

Mô hình hồi quy tuyến tính hay gọi là mô hình thống kê (model statistics) về bản chất là công thức toán biểu diễn mối quan hệ giữa các biến (các tham số trong công thức đó) với các miền giá trị của từng tham số. Chi tiết tham khảo ở link này².

Cần lưu ý là giữa hai biến \(x\) và \(y\) trong ngữ cảnh nói về sự hồi quy (regression) là ta xác định biến nào là biến độc lập, biến nào là biến phụ thuộc, còn trong ngữ cảnh khi nói về sự tương quan (correlation) thì ta chỉ xét về mức độ dao động của cả hai biến xem có cùng tăng hay cùng giảm hay không. Vì vậy phân tích tương quan là cơ sở để ta làm phân tích hồi quy.

2 Đọc kết quả của phương trình hồi quy tuyến tính

Đồ thị bên dưới cho thấy mối quan hệ giữa đường kính thân cây và thể tích thân cây có mối quan hệ tương quan tỷ lệ thuận. Công thức \(y = -36.94 + 5.07 \times x\) có \(b_0 = -36.94\) là hệ số chặn (intercept/constant) là giá trị của \(y\) khi \(x = 0\); \(b_1 = 5.07\) là hệ số góc (độ dốc/slope/coefficicent) đại diện cho mức độ tăng (hay giảm) của biến \(y\) khi \(x\) tăng lên 1 đơn vị.

Thông thường ta chỉ có thể dự đoán được biến \(y\) trong khoảng cho trước của biến \(x\) gọi là phương pháp nội suy (interpolate), các giá trị \(y\) được tính theo phương trình gọi là fitted values hay là predicted values của \(y\) (hay là giá trị \(y_{dự~đoán}\)).

Nếu cho biến \(x\) nằm ngoài khoảng giá trị của biến \(x\) ban đầu dùng xây dựng mô hình thì khi ta dự đoán biến \(y\) sẽ không chính xác (vì không đảm bảo ngoài khoảng đó thì mối quan hệ giữa \(x\) và \(y\) còn tuyến tính hay không), phương pháp ngoại suy (extrapolate) khi áp dụng cho mô hình hồi quy tuyến tính cần rất thận trọng.

There are certain guidelines for regression lines³

• Use regression lines when there is a significant correlation to predict values.
• Do not use if there is not a significant correlation.
• Stay within the range of the data. Do not extrapolate!! For example, if the data is from 10 to 60, do not predict a value for 400.
• Do not make predictions for a population based on another population’s regression line.

Khoảng chênh lệch giữa \(y_{thực~nghiệm}\) (actual observed values of \(y\)) trừ đi \(y_{dự~đoán}\) gọi là phần dư, residuals hay là sai số errors.

Phương trình đường thẳng \(y = -36.94 + 5.07 \times x\) gọi là line of best fit là một trong rất nhiều đường tuyến tính có thể có để “fit” hay là “đi qua nhiều điểm nhất với ít sai số nhất”. Việc ước lượng giá trị \(\beta_0\) và \(\beta_1\) để tìm ra giá trị \(b_0\) và \(b_1\) với p-value có ý nghĩa thống kê được gọi là quá trình xây dựng mô hình hồi quy tuyến tính.

Important

Normally we then predict values for \(y\) based on values of \(x\). This still does not mean that \(y\) is caused by \(x\).

Remember, correlation does not imply causation.

Phân tích tương quan hay phân tích hồi quy không ngụ ý là hai biến \(x\) và \(y\) có mối quan hệ nhân quả (vì có \(x\) nên mới có \(y\)) mà muốn xác định quan hệ nhân quả (cause-and-effect relationship) giữa hai biến này ta cần làm thêm một số bước phân tích về cơ chế cũng như bố trí thí nghiệm theo khối ngẫu nhiên.^4, 5

Phân tích kỹ hơn

Công thức hồi quy tuyến tính với hệ số góc \(\beta_1\) và hệ số chặn \(\beta_0\) là chỉ số của quần thể (population statistics), ký hiệu là chữ cái Hy Lạp.

\[y = \beta_0 + \beta_1 \times x\]

Điều này có nghĩa là ta có nhiều hơn 1 hệ số góc và 1 hệ số chặn để tạo thành đường hồi quy tuyến tính, hệ số góc và hệ số chặn được tìm ra bằng phương pháp bình phương tối thiểu hay bằng phương pháp nào khác được ký hiệu với dấu mũ (hat).

\[y = \hat{\beta_0} + \hat{\beta_1} \times x\]

Hay được viết theo kiểu thông thường là:

\[y = b_0 + b_1 \times x\]

\(b_0\) và \(b_1\) gọi là hệ số chặn và hệ số góc của đường best-fitting line.

3 Phân tích hồi quy tuyến tính trong R

3.1 Bước 1: Phân tích dataset

trees -> df

names(df) <- c("duong_kinh", "chieu_cao", "the_tich")

df

   duong_kinh chieu_cao the_tich
1         8.3        70     10.3
2         8.6        65     10.3
3         8.8        63     10.2
4        10.5        72     16.4
5        10.7        81     18.8
6        10.8        83     19.7
7        11.0        66     15.6
8        11.0        75     18.2
9        11.1        80     22.6
10       11.2        75     19.9
11       11.3        79     24.2
12       11.4        76     21.0
13       11.4        76     21.4
14       11.7        69     21.3
15       12.0        75     19.1
16       12.9        74     22.2
17       12.9        85     33.8
18       13.3        86     27.4
19       13.7        71     25.7
20       13.8        64     24.9
21       14.0        78     34.5
22       14.2        80     31.7
23       14.5        74     36.3
24       16.0        72     38.3
25       16.3        77     42.6
26       17.3        81     55.4
27       17.5        82     55.7
28       17.9        80     58.3
29       18.0        80     51.5
30       18.0        80     51.0
31       20.6        87     77.0

dim(df)

[1] 31  3

summary(df)

   duong_kinh      chieu_cao     the_tich    
 Min.   : 8.30   Min.   :63   Min.   :10.20  
 1st Qu.:11.05   1st Qu.:72   1st Qu.:19.40  
 Median :12.90   Median :76   Median :24.20  
 Mean   :13.25   Mean   :76   Mean   :30.17  
 3rd Qu.:15.25   3rd Qu.:80   3rd Qu.:37.30  
 Max.   :20.60   Max.   :87   Max.   :77.00  

str(df)

'data.frame':   31 obs. of  3 variables:
 $ duong_kinh: num  8.3 8.6 8.8 10.5 10.7 10.8 11 11 11.1 11.2 ...
 $ chieu_cao : num  70 65 63 72 81 83 66 75 80 75 ...
 $ the_tich  : num  10.3 10.3 10.2 16.4 18.8 19.7 15.6 18.2 22.6 19.9 ...

3.2 Bước 2: Kiểm tra giả định

Mô hình hồi quy tuyến tính cần được đáp ứng các giả định sau^{6, 7, 8, 9}

There are four principal assumptions which justify the use of linear regression models for purposes of inference or prediction:

1/ linearity and additivity of the relationship between dependent and independent variables:

1.1/ The expected value of dependent variable is a straight-line function of each independent variable, holding the others fixed.

1.2/ The slope of that line does not depend on the values of the other variables.

1.3/ The effects of different independent variables on the expected value of the dependent variable are additive.

2/ statistical independence of the errors (in particular, no correlation between consecutive errors in the case of time series data)

3/ homoscedasticity (constant variance) of the errors

3.1/ Versus time (in the case of time series data)

3.2/ Versus the predictions

3.3/ Versus any independent variable

4/ normality of the error distribution.

Việc đáp ứng được các giả định này rất quan trọng để đảm bảo kết quả ước lượng “gần” với thực tế, nếu mô hình hồi quy tuyến tính violate/không thỏa mãn các giả định trên thì kết quả ước lượng sẽ bị sai lệch (biased), dự đoán kém chính xác.

Đánh giá mô hình hồi quy tốt hay không (đáp ứng tiêu chí về độ nhạy, độ chính xác) thì ta sẽ thực hiện bước validate mô hình được đề cập ở phần sau.

4 Kiểm tra giả định về mối quan hệ giữa biến độc lập và biến phụ thuộc có là tuyến tính

4.1 Cách 1: Phân tích tương quan

Trong dataset này, ta chọn biến \(x\) là duong_kinh, biến \(y\) là the_tich, để đánh giá mức độ tuyến tính giữa hai biến này, sử dụng function cor() với method là pearson.

cor(x = df$duong_kinh,
    y = df$the_tich,
    method = "pearson")

[1] 0.9671194

4.2 Cách 2: Vẽ đồ thị scatter plot

plot(x = df$duong_kinh,
     y = df$the_tich,
     col = "red",
     pch = 1,
     xlab = "Đường kính",
     ylab = "Thể tích",
     main = "Tương quan giữa đường kính và thể tích của cây black cherry")

pairs(df)

4.3 Cách 3: Vẽ ma trận tương quan

library(corrplot)

res <- cor(df, method = "pearson",
           use = "complete.obs")

round(res, 2) -> res

corrplot(res)

5 Simple linear regression model fitting

Để thực hiện kiểm tra giả thuyết về phần dư (errors / residuals) ta cần phải xây dựng mô hình hồi quy tuyến tính thì mới có giá trị \(y_{dự~đoán}\) nhằm so sánh với \(y_{thực~nghiệm}\). Bản chất của việc đánh giá đặc điểm của phần dư này (có phân bố chuẩn không, phương sai có thay đổi hay là hằng số) chính là một phần trong các bước để validate/kiểm tra xem model này hoạt động của chính xác với giá trị thực hay không.

Ta có thể chia bộ dữ liệu ra 80% để xây dựng model (dataset train), còn lại 20% (dataset test) để kiểm tra lại xem model có dự đoán chính xác hay không. Cách tiếp cận này liên quan đến khái niệm overfitting và underfitting.

Ngắn gọn thì underfitting là model được xây dựng từ dataset train không dự đoán chính xác giá trị trong cả dataset train và dataset test, do đó model bị underfitting thì không chính xác, không sử dụng được. Còn overfitting là model dự đoán chính xác ở dataset train nhưng không dự đoán tốt ở dataset test, như vậy cũng không chính xác trong thực tế. Model có khả năng dự đoán tốt cả ở dataset train và test gọi là cần được trải qua kiểm định độ phù hợp goodness of fit.

fit <- lm(formula = the_tich ~ duong_kinh,
          data = df)

summary(fit)

Call:
lm(formula = the_tich ~ duong_kinh, data = df)

Residuals:
   Min     1Q Median     3Q    Max 
-8.065 -3.107  0.152  3.495  9.587 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -36.9435     3.3651  -10.98 7.62e-12 ***
duong_kinh    5.0659     0.2474   20.48  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.252 on 29 degrees of freedom
Multiple R-squared:  0.9353,    Adjusted R-squared:  0.9331 
F-statistic: 419.4 on 1 and 29 DF,  p-value: < 2.2e-16

par(pty = "m")

par(mar = c(5, 5, 5, 5))
par(font.axis = 2)
par(font.lab = 2)
par(cex.axis = 1.3)
par(cex.lab = 1.3)

plot(x = df$duong_kinh,
     y = df$the_tich,
     col = adjustcolor("blue", alpha.f = 0.5),
     pch = 16,
     cex = 1.5,
     xlab = "Đường kính",
     ylab = "Thể tích",
     xlim = c(0, 30),
     ylim = c(0, 100),
     xaxs = "i",
     yaxs = "i",
     las = 1)

abline(fit, 
       col = "black",
       lty = 1,
       lwd = 2)

coef(fit)[1] -> b_0
round(b_0, 2) -> b_0

coef(fit)[2] -> b_1
round(b_1, 2) -> b_1

# Tách hệ số multiple R-squared
round(unclass(summary(fit))$r.squared, 4) -> r_squared 
round(unclass(summary(fit))$adj.r.squared, 4) -> r_adj 

text(x = 10, 
     y = 80,
     labels = bquote(
       atop(
         Y == ~ .(b_0) + .(b_1) ~ "\u00D7" ~ X,
         atop(
         Multiple~R^2 == .(r_squared),
         Adjusted~R^2 == .(r_adj) 
             )
       )
                    ),
     cex = 2)

5.1 Cách tính hệ số \(Multiple~R^2\)

Hệ số R bình phương cung cấp thông tin về độ phù hợp của đường hồi quy, goodness of fit of a model¹⁰. Vì vậy đây là hệ số mang tính quyết định đến chất lượng hồi quy của model (coefficient of determination).

\(R^2 = 0.93\) có nghĩa là 93% thay đổi/biến động của \(y\) là do \(x\) gây ra khi \(x\) tăng 1 đơn vị (93% of the variation in the \(y\) values is accounted for by the \(x\) values)¹¹.

Phân biệt giữa \(Multiple~R^2\), \(R~Square\) và \(Adjusted~R^2\) (tham khảo^12, 13). Hệ số R bình phương có thể tính bằng hai cách sau:

5.1.1 Tính trực tiếp từ hệ số tương quan tuyến tính: \(Multiple~R^2 = r^2\)

\[r = \frac{\sum_{i=1}^n(x_i- \bar x)(y_i - \bar y)}{\sqrt{\sum_{i=1}^n(x_i- \bar x)^2}\sqrt{\sum_{i=1}^n(y_i- \bar y)^2}}=\frac{s_{xy}}{s_x s_y}\]

Trong đó:

• \(s_{xy}\) là hiệp phương sai (covariance) của \(x\) và \(y\)

• \(s_x\) là độ lệch chuẩn của \(x\)

• \(s_y\) là độ lệch chuẩn của \(y\)

# covariance
cov(x = df$duong_kinh,
    y = df$the_tich,
    use = "everything")

[1] 49.88812

# standard deviation of x
sd(df$duong_kinh)

[1] 3.138139

# standard deviation of x
sd(df$the_tich)

[1] 16.43785

# Pearson correlation coefficient
r_p <- cov(x = df$duong_kinh, y = df$the_tich, use = "everything") /
  (sd(df$duong_kinh) * sd(df$the_tich))

r_p

[1] 0.9671194

# Tính theo function
cor(x = df$duong_kinh,
    y = df$the_tich,
    method = "pearson")

[1] 0.9671194

# Multiple R^2
r_p^2

[1] 0.9353199

5.1.2 Tính trực tiếp từ chênh lệch giữa \(y_{thực~nghiệm}\) và \(y_{dự~đoán}\)

\[{R^2} = 1 - \frac{{sum~squared~regression~(SSR)}}{{total~sum~of~squares~(SST)}} = 1 - \frac{{\sum\limits_{i = 1}^n {({y_i} - } {{\hat y}_i}{)^2}}}{{\sum\limits_{i = 1}^n {{{({y_i} - \bar y )}^2}} }}\]

Trong đó:

• \(y_i\) là giá trị thực của \(y\) ở quan sát thứ \(i\), còn gọi là \(y_{thực~nghiệm}\)

• \(\hat y_i\) (y mũ, y hat) là giá trị dự đoán được tính từ mô hình của \(y\) ở quan sát thứ \(i\), còn gọi là \(y_{dự~đoán}\)

• \(\bar y\) (y bar) là giá trị trung bình của \(y_{thực~nghiệm}\)

Ta có 3 cách để ước lượng \(y_{dự~đoán}\) trong R:

# Tính y_dudoan từ phương trình
df$the_tich_dudoan <- coef(fit)[1] + coef(fit)[2]*df$duong_kinh 

fit$fitted.values -> dudoan_fit

predict(object = fit,
        newdata = df[, "duong_kinh", drop = FALSE]) -> dudoan_predict

# https://stackoverflow.com/questions/60357840/why-is-there-a-difference-in-linear-regression-fitted-values-and-predicted-value

names(dudoan_fit) <- NULL
names(dudoan_predict) <- NULL

identical(dudoan_fit, dudoan_predict)

[1] FALSE

all.equal(dudoan_fit, dudoan_predict)

[1] TRUE

identical(df$the_tich_dudoan, dudoan_predict)

[1] TRUE

df

   duong_kinh chieu_cao the_tich the_tich_dudoan
1         8.3        70     10.3        5.103149
2         8.6        65     10.3        6.622906
3         8.8        63     10.2        7.636077
4        10.5        72     16.4       16.248033
5        10.7        81     18.8       17.261205
6        10.8        83     19.7       17.767790
7        11.0        66     15.6       18.780962
8        11.0        75     18.2       18.780962
9        11.1        80     22.6       19.287547
10       11.2        75     19.9       19.794133
11       11.3        79     24.2       20.300718
12       11.4        76     21.0       20.807304
13       11.4        76     21.4       20.807304
14       11.7        69     21.3       22.327061
15       12.0        75     19.1       23.846818
16       12.9        74     22.2       28.406089
17       12.9        85     33.8       28.406089
18       13.3        86     27.4       30.432431
19       13.7        71     25.7       32.458774
20       13.8        64     24.9       32.965360
21       14.0        78     34.5       33.978531
22       14.2        80     31.7       34.991702
23       14.5        74     36.3       36.511459
24       16.0        72     38.3       44.110244
25       16.3        77     42.6       45.630001
26       17.3        81     55.4       50.695857
27       17.5        82     55.7       51.709028
28       17.9        80     58.3       53.735371
29       18.0        80     51.5       54.241956
30       18.0        80     51.0       54.241956
31       20.6        87     77.0       67.413183

par(pty = "m")

par(mar = c(5, 5, 5, 5))
par(font.axis = 2)
par(font.lab = 2)
par(cex.axis = 1.3)
par(cex.lab = 1.3)

plot(x = df$duong_kinh,
     y = df$the_tich,
     col = adjustcolor("blue", alpha.f = 0.5),
     pch = 16,
     cex = 1.5,
     xlab = "Đường kính",
     ylab = "Thể tích",
     xlim = c(0, 25),
     ylim = c(0, 80),
     xaxs = "i",
     yaxs = "i",
     las = 1)

abline(fit, 
       col = "black",
       lty = 1,
       lwd = 2)

coef(fit)[1] -> b_0
round(b_0, 2) -> b_0

coef(fit)[2] -> b_1
round(b_1, 2) -> b_1

# Tách hệ số multiple R-squared
round(unclass(summary(fit))$r.squared, 4) -> r_squared 
round(unclass(summary(fit))$adj.r.squared, 4) -> r_adj 

text(x = 5, 
     y = 70,
     labels = bquote(
       atop(
         Y == ~ .(b_0) + .(b_1) ~ "\u00D7" ~ X,
         atop(
         Multiple~R^2 == .(r_squared),
         Adjusted~R^2 == .(r_adj) 
             )
       )
                    ),
     cex = 2)

points(x = df$duong_kinh,
       y = df$the_tich_dudoan,
       col = adjustcolor("red", alpha.f = 0.7),
       pch = 16,
       cex = 1.5)

segments(x0 = df$duong_kinh,
         y0 = df$the_tich,
         x1 = df$duong_kinh,
         y1 = df$the_tich_dudoan,
         col = adjustcolor("darkgreen", alpha.f = 0.8),
         cex = 1,
         lwd = 3)

legend(x = "bottomright",
       y = NULL,
       legend = c("y_thực nghiệm", 
                  "y_dự đoán", 
                  "y_trung bình",
                  "phần dư (residual) dùng tính SSR", 
                  "độ lệch với trung bình dùng tính SST"),
       col = c("blue", "red", "yellow3", "darkgreen", "yellow3"),
       pch = c(16, 16, 16, NA, NA),
       lty = c(NA, NA, NA, 1, 1),
       lwd = 3,
       cex = 1.2,
       merge = FALSE)

points(x = df$duong_kinh,
       y = rep(mean(df$the_tich), dim(df)[1]),
       col = adjustcolor("yellow3", alpha.f = 0.8),
       pch = 16,
       cex = 1.5)

segments(x0 = df$duong_kinh,
         y0 = df$the_tich,
         x1 = df$duong_kinh,
         y1 = rep(mean(df$the_tich), dim(df)[1]),
         col = adjustcolor("yellow3", alpha.f = 0.5),
         cex = 1,
         lwd = 3)

Tính SSR và SST

# sum squared regression
ssr <- sum((df$the_tich - df$the_tich_dudoan)^2)

# total sum of squares
sst <- sum((df$the_tich - mean(df$the_tich))^2)

Tính \(R^2\)

1 - ssr/sst

[1] 0.9353199

identical(r_p^2, 1 - ssr/sst)

[1] FALSE

all.equal(r_p^2, 1 - ssr/sst)

[1] TRUE

# chênh lệch này do decimal tính toán của máy tính
r_p^2 - (1 - ssr/sst)

[1] -1.110223e-16

5.2 Cách tính hệ số \(Adjusted~R^2\)

Hệ số R bình phương điều chỉnh, được “adjust” theo số lượng biến \(x\) trong model, dùng để so sánh giữa các model khi có nhiều biến được đưa vào thì model nào có \(Adjusted~R^2\) cao thì dự đoán chính xác hơn¹⁴.

\[Adjusted~R^2 = 1 - \frac{{(1 - {R^2}) \times (n - 1)}}{{(n - k - 1)}}\]

Trong đó:

\(n\) là số quan sát (observation)

\(k\) là số biến (predictor variable)

Khi \(k = 1\) thì \(Multiple~R^2 = Adjusted~R^2\), khi \(k > 1\) hay là khi tăng số biến lên thì \(Multiple~R^2 > Adjusted~R^2\).

# Adjusted R^2
1 - ((1 - r_p^2) * (dim(df)[1] - 1)) / (dim(df)[1] - 1 - 1)

[1] 0.9330895

5.3 Cách tính hệ số góc \(b_1\)

The \(y\) variable is often termed the criterion variable and the \(x\) variable the predictor variable. The slope \(b_1\) is often called the regression coefficient (hệ số hồi quy) and the intercept \(b_0\) the regression constant (hằng số hồi quy)¹⁵.

\[slope = \frac{{rise}}{{run}} = \frac{{dy}}{{dx}} = \frac{{\Delta y}}{{\Delta x}} = \frac{{{y_2} - {y_1}}}{{{x_2} - {x_1}}}\]

Tham khảo cách tính thủ công ở đây^16, 17

5.3.1 Tính từ hệ số tương quan và độ lệch chuẩn

\[b_1 = r \times \left( {\frac{{{s_y}}}{{s{}_x}}} \right)\]

Trong đó:

• \(b_1\) là hệ số góc (slope)

• \(r\) là hệ số tương quan tuyến tính theo Pearson

• \(s_x\) là độ lệch chuẩn của \(x\)

• \(s_y\) là độ lệch chuẩn của \(y\)

b_1a <- r_p * (sd(df$the_tich) / sd(df$duong_kinh))
b_1a

[1] 5.065856

5.3.2 Tính từ giá trị \(x\) và \(y\)

\[{b_1} = \frac{{n \times (\sum {xy) - (\sum {x)} } \times (\sum {y)} }}{{n \times (\sum {{x^2}) - {{(\sum x )}^2}} }}\]

b_1b <- (dim(df)[1] * sum(df$duong_kinh*df$the_tich) - sum(df$duong_kinh) * sum(df$the_tich)) /
  (dim(df)[1] * sum(df$duong_kinh^2) - sum(df$duong_kinh)^2)
b_1b

[1] 5.065856

So sánh với hệ số góc tính bằng R

coef(fit)[2] -> b_1
b_1

duong_kinh 
  5.065856 

5.4 Cách tính hệ số chặn \(b_0\)

5.4.1 Tính từ hệ số góc và trung bình của \(x\) và \(y\)

\[b_0 = \bar{y} - b_1 \times \bar{x}\]

Trong đó:

• \(b_0\) là hệ số chặn (intercept)

• \(b_1\) là hệ số góc (slope)

• \(\bar{y}\) là trung bình của các giá trị \(y_{thực~nghiệm}\)

• \(\bar{x}\) là trung bình của các giá trị \(x_{thực~nghiệm}\)

b_0a <- mean(df$the_tich) - (b_1a * mean(df$duong_kinh))
b_0a

[1] -36.94346

5.4.2 Tính từ giá trị \(x\) và \(y\)

\[{b_0} = \frac{{(\sum y) \times (\sum {x^2}) - (\sum x) \times (\sum xy)}}{{n \times (\sum {x^2}) - {{(\sum x )}^2}}}\]

b_0b <- ((sum(df$the_tich) * sum(df$duong_kinh^2)) - (sum(df$duong_kinh) * sum(df$duong_kinh * df$the_tich))) /
  (dim(df)[1] * sum(df$duong_kinh^2) - sum(df$duong_kinh)^2)
b_0b

[1] -36.94346

So sánh với hệ số chặn tính bằng R

coef(fit)[1] -> b_0
b_0

(Intercept) 
  -36.94346 

5.5 Cách tính độ lệch chuẩn của phần dư

5.5.1 Tính trực tiếp từ chênh lệch giữa \(y_{thực~nghiệm}\) và \(y_{dự~đoán}\)

\[e = errors = residuals = y_{thực~nghiệm} - y_{dự~đoán}\]

df$residuals <- df$the_tich - df$the_tich_dudoan
df

   duong_kinh chieu_cao the_tich the_tich_dudoan  residuals
1         8.3        70     10.3        5.103149  5.1968508
2         8.6        65     10.3        6.622906  3.6770939
3         8.8        63     10.2        7.636077  2.5639226
4        10.5        72     16.4       16.248033  0.1519667
5        10.7        81     18.8       17.261205  1.5387954
6        10.8        83     19.7       17.767790  1.9322098
7        11.0        66     15.6       18.780962 -3.1809615
8        11.0        75     18.2       18.780962 -0.5809615
9        11.1        80     22.6       19.287547  3.3124528
10       11.2        75     19.9       19.794133  0.1058672
11       11.3        79     24.2       20.300718  3.8992815
12       11.4        76     21.0       20.807304  0.1926959
13       11.4        76     21.4       20.807304  0.5926959
14       11.7        69     21.3       22.327061 -1.0270610
15       12.0        75     19.1       23.846818 -4.7468179
16       12.9        74     22.2       28.406089 -6.2060887
17       12.9        85     33.8       28.406089  5.3939113
18       13.3        86     27.4       30.432431 -3.0324313
19       13.7        71     25.7       32.458774 -6.7587739
20       13.8        64     24.9       32.965360 -8.0653595
21       14.0        78     34.5       33.978531  0.5214692
22       14.2        80     31.7       34.991702 -3.2917021
23       14.5        74     36.3       36.511459 -0.2114590
24       16.0        72     38.3       44.110244 -5.8102436
25       16.3        77     42.6       45.630001 -3.0300006
26       17.3        81     55.4       50.695857  4.7041430
27       17.5        82     55.7       51.709028  3.9909717
28       17.9        80     58.3       53.735371  4.5646292
29       18.0        80     51.5       54.241956 -2.7419565
30       18.0        80     51.0       54.241956 -3.2419565
31       20.6        87     77.0       67.413183  9.5868168

df$residuals

 [1]  5.1968508  3.6770939  2.5639226  0.1519667  1.5387954  1.9322098 -3.1809615 -0.5809615  3.3124528  0.1058672
[11]  3.8992815  0.1926959  0.5926959 -1.0270610 -4.7468179 -6.2060887  5.3939113 -3.0324313 -6.7587739 -8.0653595
[21]  0.5214692 -3.2917021 -0.2114590 -5.8102436 -3.0300006  4.7041430  3.9909717  4.5646292 -2.7419565 -3.2419565
[31]  9.5868168

unclass(summary(fit))$residuals

         1          2          3          4          5          6          7          8          9         10 
 5.1968508  3.6770939  2.5639226  0.1519667  1.5387954  1.9322098 -3.1809615 -0.5809615  3.3124528  0.1058672 
        11         12         13         14         15         16         17         18         19         20 
 3.8992815  0.1926959  0.5926959 -1.0270610 -4.7468179 -6.2060887  5.3939113 -3.0324313 -6.7587739 -8.0653595 
        21         22         23         24         25         26         27         28         29         30 
 0.5214692 -3.2917021 -0.2114590 -5.8102436 -3.0300006  4.7041430  3.9909717  4.5646292 -2.7419565 -3.2419565 
        31 
 9.5868168 

Phương sai của phần dư

\[s_{residuals}^2 = \frac{{\sum\limits_{i = 1}^n {e_i^2} }}{{n - 2}}\]

Trong đó:

• \(e\) là phần dư, chênh lệch giữa \(y_{thực~nghiệm} - y_{dự~đoán}\) ở từng quan sát

var_error <- sum(df$residuals^2) / (dim(df)[1] - 2) # 2 ở đây là hai biến x và y
var_error

[1] 18.0794

Độ lệch chuẩn của phần dư¹⁸ hay còn gọi là sai số chuẩn ước lượng standard error of estimate hay là residual standard error.

sd_error <- sqrt(var_error)
sd_error

[1] 4.251988

5.6 Tính theo công thức

\[{s_{residuals}} = {s_y} \times \sqrt {1 - {r^2}} \times \sqrt {\frac{{n - 1}}{{n - 2}}}\]

sd_error_ok <- sd(df$the_tich) * sqrt(1 - r_p^2) * sqrt((dim(df)[1] - 1) / (dim(df)[1] - 2))
sd_error_ok

[1] 4.251988

Theo công thức này ta thấy rõ ảnh hưởng của số lượng mẫu quan sát đến sai số khi dự đoán. In this case the error is less than 2% when \(n > 26\) and less than 1% when \(n > 51\). However, it is 10% or larger when \(n < 8\)

So sánh với cách tính residual standard error bằng R

# Residual standard error: 4.252 on 29 degrees of freedom
unclass(summary(fit))$sigma

[1] 4.251988

Nghĩa là kết quả dự đoán có độ lệch chuẩn so với giá trị thực nghiệm là 4.25

Tip

The standard error is small when the correlation is high. This increases the accuracy of prediction.

When we consider multiple distributions it is often assumed that their standard deviations are equal.

This property is called homoscedasticity.

We often consider the conditional distribution or distribution of all y scores with the same value of x.

If we assume these conditional distributions are all normal and homoscedastic, we can make probabilistic statements about the predicted scores.

The standard deviation we use is the standard error calculated above.

6 Kiểm tra giả định về phần dư của mô hình tuyến tính

6.1 Phần dư có phân bố chuẩn?

plot(fitted(fit), resid(fit))

https://sphweb.bumc.bu.edu/otlt/MPH-Modules/BS/R/R5_Correlation-Regression/R5_Correlation-Regression5.html

có biện luận kết quả

7 Tài liệu tham khảo

7.1 Regression

1. https://www.personality-project.org/r/r.commands.html

2. https://www.statology.org/linear-regression-assumptions/

3. https://online.stat.psu.edu/stat500/lesson/9/9.2/9.2.3

4. https://analystprep.com/study-notes/cfa-level-2/assumptions-of-the-simple-linear-regression-model/

5. https://datatab.net/tutorial/pearson-correlation

6. https://www.geo.fu-berlin.de/en/v/soga-r/Basics-of-statistics/Descriptive-Statistics/Measures-of-Relation-Between-Variables/Correlation/index.html

7. https://stackoverflow.com/questions/12201439/is-there-a-difference-between-the-r-functions-fitted-and-predict

7.2 R

1. https://stackoverflow.com/questions/19053440/r-legend-with-points-and-lines-being-different-colors-for-the-same-legend-item

Footnotes

1. https://en.wikipedia.org/wiki/Ordinary_least_squares

2. https://www.stat.uchicago.edu/~pmcc/pubs/AOS023.pdf

3. https://www.andrews.edu/~calkins/math/edrm611/edrm06.htm

4. https://pedermisager.org/blog/why_does_correlation_not_equal_causation/

5. https://en.wikipedia.org/wiki/Bradford_Hill_criteria

6. https://stataguide.wordpress.com/2020/04/19/cac-gia-dinh-cua-mo-hinh-hoi-quy-tuyen-tinh/

7. https://sphweb.bumc.bu.edu/otlt/MPH-Modules/BS/R/R5_Correlation-Regression/R5_Correlation-Regression4.html

8. https://www.ncl.ac.uk/webtemplate/ask-assets/external/maths-resources/statistics/regression-and-correlation/assumptions-of-regression-analysis.html

9. https://people.duke.edu/~rnau/testing.htm

10. https://www.ncl.ac.uk/webtemplate/ask-assets/external/maths-resources/statistics/regression-and-correlation/coefficient-of-determination-r-squared.html

11. https://www.ncl.ac.uk/webtemplate/ask-assets/external/maths-resources/statistics/regression-and-correlation/coefficient-of-determination-r-squared.html

12. https://maths.uel.edu.vn/Resources/Docs/SubDomain/maths/TaiLieuHocTap/ToanUngDung/h_s_r_bnh_phng_v_r_bnh_phng_hiu_chnh.html

13. https://www.statology.org/multiple-r-vs-r-squared/

14. https://www.statology.org/multiple-r-vs-r-squared/

15. https://www.andrews.edu/~calkins/math/edrm611/edrm06.htm

16. https://www.statisticshowto.com/probability-and-statistics/regression-analysis/find-a-linear-regression-equation/

17. https://www.dummies.com/article/academics-the-arts/math/statistics/how-to-calculate-a-regression-line-169795/

18. https://www.andrews.edu/~calkins/math/edrm611/edrm06.htm